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ABSTRACT 

This  paper  discusses  several  approaches  to  handling  the 
filtering  and  smoothing  of  time  related  data.  The  specific 
case  addressed  is  the  Antitank  Missile  Test  (ATMT)  conducted 
at  US  Army  Combat  Developments  Experimentation  Conmand  using 
the  Range  Measuring  System.  Methods  discussed  include  difference 
quotients,  Kalman  filters,  spline  functions,  least  squares 
approximation,  convolution  filters,  and  minimim  maximun  error 
techniques. 
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FILTERING  AND  SMOOTHING  OF  TIME  RELATFO  DATA 


! 

by  David  W.  Bash,  PhD 

I.  INTRODUCTION. 

, a.  This  paper  discusses  several  approaches  to  handling  data  recorded 

by  the  Range  Measuring  System  (RMS)  at  US  Army  Combat  Developments 
Experimentation  Command  (CDEC)  to  obtain  velocity  and  acceleration 
estimates  for  vehicles.  The  discussion  also  applies  to  any  data  when 
approximate  X,  Y,  Z coordinates  as  a function  of  time  are  given. 

b.  fs  in  most  data  handling  situations,  the  computation  of  velocity 
and  acceleration  is  highly  dependent  upon  the  physical  setting,  the 

i accuracy  of  the  measurements,  and  the  desired  application  of  the  results. 

This  paper  was  prepared  as  a result  of  investigations  made  during  analysis 
! of  the  Antitank  Missile  Test  (ATMT)  conducted  at  CDEC  and  is  illustrated  with 

data  from  that  study.  The  ATiyTT  analysis  was  interested  in  following  the 
motions  of  the  vehicle  closely  in  order  to  ascertain  the  effects  of 
evasive  maneuvers  of  targets  on  certain  weapon  systems. 

c.  The  accurate  position,  velocity,  and  acceleration  of  a vehicle 
may  be  required  for  many  test  purposes.  This  paper  shows  that  some 
approaches  to  calculating  velocity  and  acceleration  yield  erroneous 

j results.  ATMT  used  RMS  data  in  terms  of  time,  X-coordinate,  Y-coordinate, 

and  Z-coordinate  for  a variety  of  target  vehicles  traversing  several 
general  paths  on  different  terrains  and  following  various  operating 
instructions.  RMS  data  give  the  X-,  Y-,  and  Z-coordinates  in  meters 
approximately  10  times  per  second.  Usually,  but  not  always,  the  time 
differences  between  consecutive  measurements  in  ATMT  were  0.1  + 0.01 
i seconds.  The  error  in  coordinates  was  reported  to  have  a standard 

I deviation  of  1 meter,  a mean  of  zero,  and  a normal  distribution.  How- 

i ever,  the  data,  as  received,  had  already  been  processed  by  a Kalman  filter 

i so  that  actually  a maximum  likelihood  estimate  of  the  true  location  was 

being  reported.  Thus,  the  three  major  errors  in  computing  the  velocity 
; and  acceleration  of  the  vehicle  were  the  following: 

(1)  The  data  error  in  location,  referred  to  in  this  paper  as 
the  inherent  error. 

(2)  The  error  in  computer  computation,  referred  to  as  the 
roundoff  error. 

(3)  The  error  by  using  a finite  approximation  formula,  referred 
to  as  the  truncation  error. 


i 
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d.  Much  of  the  discussion  in  the  following  paragraphs  is  with  respect 
to  just  one  coordinate.  However,  all  coordinates  need  to  be  considered, 
and  the  magnitude  and  direction  of  velocity  and  acceleration  vectors 
found,  since  maneuvers  may  impart  more  velocity  and  acceleration  than 
the  vehicle  can  perform  in  a straight  line.  This  paper  discusses  the 
strengths  and  weaknesses  of  several  methods  for  application.  The 
references  listed  at  the  end  of  the  paper  give  more  complete  Implementation 
information;  the  algorithms  included  here  are  intended  merely  to  Illustrate 
the  various  methods. 

2.  DIFFERENCE  QUOTIENT  APPROACH. 

a.  Standard  Approach.  The  first  technique  discussed  is  the  straight- 
forward approach  of  using  the  definition  of  first  and  second  derivatives 
to  conipute  difference  quotients.  Table  1 presents  a section  of  the 
results  of  approximating  the  first  derivative  of  x(t^)  by: 


x'(t^) 


and  the  second  derivative  by: 


x(tj  + j)  - ^(tj) 
^i  . 1 ■ S 


(Eq  1) 


x’(t.  + ,)  - x'(tJ 
;'-(tJ  = 2- 


(Eq  2) 


^i  + 2 ■ ^i 


Table  1.  Standard  difference  quotients 


t^:  (sec) 

x(t^):  (meters) 

x'(tj: 

(meteri/sec) 

x"(t.): 

(meters/sec‘^) 

37010.76 

55195.9 

5.00 

25.25 

37010.82 

55196.2 

7.78 

-44.44 

37010.91 

55196.9 

4.44 

17.28 

37011.00 

55197.3 

6.00 

- 5.74 

37011.10 

55197.9 

5.45 



0.96 

Accolorations  larger  than  straight-line  speedup  or  braking  could  occur 
in  a sharp  turn,  but  values  of  at  least  4 g (one  q is  the  acceleration 
duc'  to  (jravity  of  about  0.0  meters  per  second’)  are  unreasonable  for 
the  vehicles  being  considered,  as  are  changes  in  acceleration  of  more 
than  7 g in  one- tenth  of  a second.  The  main  source  of  error  in  this 
approach  is  the  inherent  error.  A small  error;  e.g.,  one-half  standard 
deviation  or  one-half  meter  in  just  one  x(t)  value  of  the  first  two 
entries  in  table  1,  could  cause  100  percent  error  in  the  x'(t)  value. 
Actually,  the  inherent  error  in  the  acceleration  is  proportional  to  the 
error  in  the  data  divided  by  the  time  difference  squared,  roughly  a 
factor  of  100  times  the  data  error. 

b.  Point  Skipping  Approach.  One  approach  to  controlling  the  error 
difference  quotients  would  be  to  increase  the  time  difference  between 
consecutive  points  (i.e.,  skip  some  data  points)  to  hold  down  the  inherent 
error.  A method  like  this  could  be  sound,  at  least  heuristically.  High 
frequency  random  measuring  errors  may  be  prevented  from  affecting  the 
differences  used  in  computing  difference  quotients.  However,  in  an 
application  interested  in  evasive  maneuvers,  no  predetermined,  large- 
size  skipping  rule  can  guarantee  that  evasive  maneuvers  are  not  missed. 
Table  2 displays  the  difference  quotients  for  a small-size  skipping  rule. 
Again,  these  values  are  too  rough  to  use  in  most  applications.  Another 
source  of  error  is  that,  since  the  velocities  and  accelerations  are  not 
constant  functions  of  time,  the  truncation  error  is  proportional  to  a 
power  of  the  time  difference  (see  chapter  4 of  Ralston,  reference  8). 

Thus,  increasing  the  time  differences  may  increase  the  truncation  error. 

Table  2.  Skipping  points  difference  quotients 


t^(sec) 

x(t^)  (meters) 

*'(tj 

(me tens /sec) 

x"(t.) 

(meters/dec*^) 

37010.76 

55195.9 

5.78 

- 1.93 

37011.21 

55198.5 

4.81 

- 4.68 

37011.75 

55201.1 

2.50 

0.00 

37012.23 

55202.2 

2.50 

- 3.50 

37012.79 

55203.7 

0.68 

0.82 
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c.  Int^polation.  Another  approach  is  to  Interpolate  the  data  with 
some  type  of  function  and  compute  the  first  and  second  derivatives  of 
the  function.  Usually,  a polynomial  is  tried  (e.g.,  Lagrange  polynomials, 
Uermite  polynomials),  but  rational  functions  or  transcendental  functions 
[K)ssibly  could  be  used.  One  trouble  with  this  approach  is  in  choosing 
the  degree  of  other  parameters  in  the  expression.  If  the  degree  is  low 
(i.e.,  only  a few  points  at  a time  are  interpolated),  the  polynomials  will 
follow  the  errors  of  the  data  and  the  derivatives  of  the  polynomials  will 
be  similar  to  the  difference  quotients  discussed  previously.  If  the 
degree  is  high  (i.e.,  many  points  at  a time  are  interpolated)  the  poly- 
nomials will  more  than  likely  have  extreme  and  unreasonable  variations 
(see  page  354  of  Hamning,  reference  6).  Again,  one  might  try  employing 
only  selected  data  points,  but  the  same  problems  occur  here  as  in  the 
point  skipping  approach.  Thus,  these  direct  approaches  almost  always 
yield  questionable  results. 

3.  KALfiAN  FILTER.  Subparagraphs  3a  and  3b,  below,  discuss  the  use  of 
the  Kalman  filter  in  filtering;  i.e.,  changing  the  data  point  by  point 
as  they  are  presently  being  received  or  considered.  Subparagraph  3c 
discusses  smoothing;  i.e.,  looking  at  and  changing  past  values. 

a.  Filtering.  If  the  errors  are  distributed  normally  with  mean 
zero,  the  Kalman  filter,  which  is  a weighted-least-squares  estimate  of 
the  true  coordinates,  gives  a maximun-l ike li hood  or  minimun-variance 
estimate  of  the  true  location.  Matrix  vector  notation,  with  a super- 
script T for  the  transpose,  is  used  in  the  following  discussion.  Let: 

x^  = actual  position  vector 

£ = measured  position  vector 

i = estimated  position  vector 

H = a 3x3  matrix  stating  the  assiwed  relation  of  actual  to  measured 
V = random  error  vector 


z = Hx  + V (Eq  3) 

r(y) 

M 


= a 3x3  variance-covariance  matrix  for  v = E((£  - H£)(^  - Hx)^) 

= a 3x3  variance-covariance  matrix  for  >1  - ^ = E((x  - x)(x  - 2)^) 
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The  objective  in  Kalman  filtering  is  to  find  as  the  value  of  j<  that 
minimizes  the  quadratic  form: 


J = j^(x  - M'*  (x  - x)  + (/  - Hx_)^  R'*  (z  - Hx)j  (Eq  'i) 

The  X is  found  by  standard  differential  argunents.  Note  that  errors  x - x 
and  z - are  weighted  by  the  prior  expected  values  of  {x  - >^{x  - 
and  Jz  - Hx)(£  - Hx)^,  respectively. 

b.  Multistage  Filtering.  The  situation  is  more  complicated  in  multi- 
stage processes,  as  exemplified  by  the  ATMT  study.  In  the  following 
discussion,  the  notation  is  changed  slightly  to  acconmodate  a linear  first 
order  prediction  step  for  one  model  of  RMS  data.  Assuming  the  best  least- 
squares  estimate  for  the  (i  - l)th  point,  the  best  least-squares  estimate 
for  the  ith  point  is  stated.  Let: 


Xi  = best  least-squares  estimate  for  position-velocity  vector 
= (Xj.  y,,  X,,  z,) 


♦ ^ j = transition  matrix  from  state  of  time  t.  _ j to  time  t. 


, a 3x3  identity  matrix 


, a 6x6  matrix 
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= random  error  vector 


(Eq  5) 


= measured  position  vector  * (x^  , z^ , 0,  0,  0)  ■ •*  ^ 

= expected  value  matrix  for  (£  - - ^1-i^^  " ^ili^ 


^i  = * i - 


M.  - 4.  , 


1 X i _ 1 = predicted  position-velocity  vector 
T 


^ - l^i 
Pi  = (M.  ■ W H.\  ■ ^ H.) 


i - 1 i - 1 i - 1 

T„  - 1 ..  X - 1 


(Eq  6) 


r.  xT 


= expected  value  matrix  for  (x  - ^)(x  - 


Kj  = P^h/r,  - ' 


(Eq  7) 


Then  the  best  least-squares  estimate  is  a standard  differential  change 
predicted  value  x-  plus  a correction  term  proportional  to  the  difference 
between  measured  And  predicted  value: 


ii  ' ii  * (if  - i i 1. 


(Eq  8) 


For  computational  purposes,  it  is  helpful  if  the  dt.'s  are  small  compared 
to  1.  Then  the  asstmption  is  made  that  dt-  = 0,  which  implies  that  , 
is  the  identity  in  equation  6.  Also,  it  is  assuned  that  the  errors  ' " ^ 
in  the  coordinates  are  independently  and  identically  distributed  so  that 
and  are  all  diagonal.  Then  i.,  with  these  assumptions,  is  only 
approximately  the  maximum  likelihood  estimate.  However,  these  assumptions 
may  be  needed  to  speed  up  calculations  to  run  in  real  time,  as  was  desired 
in  the  original  filtering  of  the  ATMT  data  (reference  3).  Chapter  12  of 
Bryson  and  Ho  (reference  4)  contains  the  details  of  the  above  derivation, 
and  the  BDM  docunentation  (reference  3)  includes  the  practical  simplifications. 
Reference  2 describes  a similar  model  accommodating  acceleration,  and 
chapter  4 of  the  text  by  McGarty  (reference  7)  gives  another  readable 
account.  Table  3 is  a partial  output  from  a proqramned  version  of  the 
algorithm  in  the  McGarty  text. 
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Table  3.  Kalman  filter  example 


t,:  (sec) 

x(t,):  (meters) 

1 

x"(t.)* 

(meters/iec^) 

34130.13 

56464.9 

6.06 

1.47 

34130.24 

56465.5 

6.22 

1.48 

34130.33 

56466.2 

6.35 

1.47 

34130.43 

56466.7 

6.45 

1.45 

34130.59 

56467.4 

6.61 

1.46 

c.  Smoothing.  After  all  the  data  have  been  recorded  and  filtered, 
then  one  may  start  with  the  last  value  Xu  and  work  forward,  smoothing  the 
data  by  minimizing  a similar  but  more  explicated  functional.  The  results 
that  follow,  and  some  necessary  initializing  conditions,  are  derived  in 
chapter  13  of  Bryson  and  Ho.  Let: 


^.i/N  " least-squares  estimate  to  i 

N points. 


^th 


point  given 


(Eg  9) 


Then: 


iii/N  ^ ii  ■ ^i 


(i  i + 


1 ■ - i 


+ 1/N^ 


(Eq  10) 


where: 


Pi*, 


- 1 

M + 1* 


(Eq  11) 


The  objective  in  Kalman  filtering  and  smoothing  is  to  change  the  dependent 
variable  values  to  be  nice  in  a statistical  sense.  Knowledge  of  the 
variance-covariance  matrices  discussed  previously  is  necessary  before 
doing  Kalman  filtering  and/or  smoothing.  Often,  the  variances  are  adjusted 
to  yield  nice  output. 
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4.  POLYNOMIAL  SPLINE  FUNCTIONS. 
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a.  Standard  Approach.  Spline  function  approximation  is  the  inter- 
polation  of  data  by  polynomials  of  a given  degree,  each  polynomial  defined 
on  distinct  intervals  between  consecutive  data  points.  These  polynomials 
are  determined  so  that  the  resulting  function  defined  on  the  entire  domain 
has  continuous  derivatives  of  all  orders  up  to  but  excluding  the  degree 

of  the  polynomial  (even  at  the  data  points).  A property  of  interest  in 
spline  approximation  is  that,  using  polynomials  of  degree  2k-l,  the 
integral : 

b 

i (f^'‘^(t))2dt  (Eg  12) 

a 

is  a miniimw  over  all  functions  with  the  same  interpolation  and  smoothness 
property.  In  the  case  of  cubic  polynomials  (the  most  often  used  degree) 
the  second  derivative,  or  accelerations,  is  minimized  in  this  sense.  Tltis 
property  of  minimal  acceleration  may  not  be  too  severe  even  if  the  project 
is  mainly  interested  in  actual  accelerations.  The  reasons  include  the 
fact  that  the  function  must  have  sufficiently  large  second  derivative  in 
order  to  interpolate  the  data  points  and  most  physical  systems  operate  in 
a conservation  of  work  (proportional  to  acceleration)  principle.  The 
spline  approach  can  be  advantageous  if  the  values  of  the  first  and  second 
derivatives  are  needed  at  points  other  than  data  points.  In  the  derivation 
of  these  polynomials  nice  computational  formulas  are  developed  for  the 
derivatives  at  any  domain  value  (subparagraph  4b,  below).  This  capability 
might  be  needed  if  the  data  points  are  rather  sparse  compared  to  RMS  data 
and  the  unreasonable  variations  mentioned  in  paragraph  2 are  to  be  avoided. 
However,  in  order  to  obtain  some  of  the  parameters  ir  the  formulas,  it  is 
necessary  to  invert  an  n by  n matrix  where  n is  the  number  of  data  points 
being  fitted.  The  matrix  is  fairly  diagonal  dominant,  but  care  must  still 
be  used  to  keep  the  roundoff  error  under  control  if  the  number  of  points 
is  large.  These  problems  can  be  alleviated  if  the  fitting  is  done  in 
segments  or  if  only  selected  points  are  used;  however,  the  first  approach 
may  destroy  some  desired  continuity  and  the  second  may  miss  some  evasive 
maneuvers  as  discussed  in  subparagraph  2c.  Also,  the  inherent  error  may 
become  serious.  The  algorithms  deriving  the  parameters  essentially  compute 
first  and  second  order  difference  quotients  and  thus  are  open  to  the 
criticisms  raised  in  subparagraphs  2a  and  2b. 

b.  Formulas  for  Cubic  Spline  Interpolation.  Given  the  following: 
t^ , 1 “"I  , . . . . , n 


P li—li  •••  • fH 


the  following  assumptions  are  made: 


x(t)  = x^(t),  £ t ^ t^  _ i«l,  ....  n - 1 

where  the  x^.(t)  are  third  degree  polynomials  in  t 


XjCtj) 

(Eq  14) 

x,(t,) 

X 

II 

X 

11 

_ ^ 1 j . . . • n ~ 1 

l(‘n> 

c 

X 

II 

1 

I 

1 

.(t,) 

■ 

i=2,  . . . , n - 1 

l(t,) 

= ^(t,) 

i=2,  ....  n - 1 

The  functional  values,  i=l,  n - 1,  are: 


\ 


S i " ♦ I 


■'r.  - ' "n 


x_-x, 


«n-l(V= 


•n  - 1 


’■-.  - l<‘n)  = '‘n 


where  the  c^. , i=l,  . . . , n are  the  solutions  to: 


AC  = D 


(Eq  16) 


The  entries  in  the  matrix  A and  vector  D are: 


i»J“l»  • • • »n 


D = (d.) 


a 
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J and  K are  approximations  to  x{t)  at  tj  and  t^,  respectively.  Further 
information  may  be  obtained  in  chapter  20  of  Hamming.  Figures  1,  2,  and 
3 show  the  plots  of  position,  velocity,  and  acceleration,  respectively, 
using  the  cubic  spline  approach  to  smoothing  using  every  eleventh  data 
point.  Note  the  acceleration  is  piecewise  linear  and  continuous  as  it 
should  be  but  still  is  realistic. 

.c.  Spline  Approximation  with  Least  Squares  Approach.  This  approach  is 
designed  1b  minimize  the  inherent  error  problem  mentioned  in  the  previous 
paragraph.  The  usual  approach  minimizes  the  weighted  sum  of  a smoothness 
term  and  the  sum  of  the  squares  of  the  deviations  between  observed  values 
and  smoothed  values.  Let: 


E(x)  = D (x(t.)  - x.)^ 
i=l  ' ’ 


(Eq  17) 
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Figure  1.  Position  plot 


S(x)  = 


(x^t))^  dt 


(Eq  18) 


b 

I 

a 


n = 2k  - 1,  k integer 


(Eq  19) 


0 < p < 1 


Then  the  smoothing  spline  function  x(t)  for  the  data  points  (t.,  x^),  using 
polynomials  of  degree  n,  is  the  function  that  minimizes: 


S(x)  + pE(x) 


(Eq  20) 


The  factor  p weights  the  relative  importance  of  the  two  terms.  Given  p, 
to  solve  for  x(t)  it  is  required  to  solve  a system  of  equations  similar  to 
the  systems  in  paragraphs  4a  and  5c.  Hence,  the  cautions  with  respect  to 
computation  apply  here.  However,  there  is  no  interpolation  now,  but  the 
first  and  second  derivations  will  be  smoother.  Some  of  the  reasons  that 
this  method  may  not  be  desirable  as  applied  to  RMS  data  are  the  fact  that 
derivatives  at  non-data  points  are  not  needed  and  there  are  too  many  points 
to  use  all  at  one  time,  requiring  doing  the  analysis  in  segments  (destroying 
some  continuity)  or  using  selected  points  (raising  concern  about  missing  a 
maneuver).  More  detail  on  implementation  may  be  found  in  Schimaker, 
reference  10. 

5.  LEAST  SQUARES  APPROXIMATION  (FILTERS). 

a.  Equal  Spacing.  If  the  independent  variable,  time,  is  equally  spaced 
(i.e.,  the  difference  between  consecutive  times  is  constant),  then  nice 
computational  fonnulas  can  be  used  to  estimate  first  and  second  derivatives. 
For  example,  if  x(t)  is  a position  function  of  time,  then  y'(t.  ) is 
approximately:  ^ 

- -j7  |^x(t|^  + 2^  ■ ^^^k  - 2O  * [^^^k  + 1^  ■ ’^^^k  - 1^] 

(Eq  21) 

where  the  t.  are  equally  spaced.  This  expression  is  exact  for  all  poly- 
nomials of  degree  2 or  less  if  the  differences  of  consecutive  independent 
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vdri.ihlos  arc  one.  The  restriction  of  differences  all  one  can  be  overcome 
by  .i|)|iro|)r iate  scaling  in  the  independent  variable.  This  type  of  formula 
can  be  derived  using  either  Fourier  series  or  differentiating  a least-squares 
best  fit  (wlynomlal,  a special  case  of  subparagraphs  b,  c,  and  d,  below. 

The  polynomial  Is  the  polynomial  of  degree  n (where  n + 1 is  less  than 
the  number  of  data  points  being  fit),  which  minimizes  the  sum  of  the 
squares  of  the  difference  in  the  observed  value  and  polynomial  value  at 
the  points  being  fitted. 

b.  Unequ^al  Spacing.  This  was  the  situation  encountered  in  the  ATMT 
study.  The  formuUs  in  this  case  are  not  as  simple  as  those  in  subparagraph 
a,  above.  The  idea  again  is  to  fit  a polynomial  of  degree  n: 

x(t)  = ag  + aj  t + . . . . + a^  _ jt"  ■ ^ + a^t"  (Eq  22) 


to  m points,  m > n + 1,  which  minimizes  the  sum  of  the  squares  of  the 
differences  between  the  observed  value  and  x(t)  at  the  m points.  The  n 
and  m are  determined  from  considering  the  application.  For  example,  the 
degree  of  polynomial  n may  be  heuristically  set  by  matching  the  nunber 
of  major  acceleration  changes  during  the  interval  of  m points  with  the 
number  of  sign  changes  possible  with  the  second  derivative  of  an  nth  degree 
polynomial  on  the  interval  of  m points.  For  a given  set  of  data  and  n, 
the  larger  the  m (the  number  of  points)  the  "smoother"  the  polynomial 
values  and  derivatives  will  be.  Thus,  m must  be  large  enough  to  smooth 
the  values  so  that  the  velocities  and  accelerations  are  within  the 
capability  of  the  vehicle  system.  An  upper  bound  on  m possibly  can  be 
determined  by  checking  if  values  are  smoothed  below  known  actual  performance 
or  by  an  autocorrelation  program  on  the  data,  with  trend  and  any  discernible 
cycles  removed  to  determine  when  the  errors  become  essentially  linearly 
independent.  In  the  ATMT  analysis,  the  degree  used  was  3 and  the  nunber 
of  points  was  between  9 and  51.  Different  values  would  probably  be  needed 
for  aircraft  or  other  vehicles.  The  general  approach  with  a large  quantity 
of  data  is  to  fit  a polynomial  to  each  t.  , where  (m  + l)/2  < k £ maximum 
number  of  points  minus  (m  - l)/2,  and  use  the  coefficients  to  compute  the 
smoothed  position,  velocity,  and  acceleration  in  each  component  as  a 
function  of  time.  Thus,  the  problem  in  the  final  determination  of  m is 
to  balance  realistic  velocity  and  acceleration  with  an  analysis  of  the 
distances  between  the  observed  and  smoothed  positions.  A complete 
statistical  analysis  of  residuals  would  require  an  immense  amount  of 
computer  time. 


c.  IX‘r iv.it  ioti  op  Co<*>  I icierits.  Ihc  followinfi  nol/it  ion  is  usimI  in 
this  discussion: 

x(t)  = the  actual  x-coordinate  (unknown) 
x(t)  = the  observed  x-coordinate 
x(t)  = the  smoothed  x-coordinate 

x(t)  + . e a^t^  (E, 

m = number  of  points  to  be  fitted 
L = (m  - l)/2 

t.  = the  time  in  question:  (m  + l)/2  £ k ^ maximun  number  of 
points  minus  (m  - l)/2 


(Eq  23) 


x(t.)  = z a.(t.)^ 
j=0  J 


1 1 

X component  of  velocity  at  t.  = Z a.j(t. ) 


j - 1 


X component  of  acceleration  at  t.  = Z a.j(j  - l)(t.) 

K j^2  J 


.]•  - 2 


The  objective  in  least  squares  fit  is  to  determine  the  coefficients  a 
that  minimize:  P 

k + L r n .T  „ 

H(a  , a,,  . . . a ) = E x(t.)  - Z a.(t.)J  ^ (Eq  24) 

u ^ " i=k  - L L ^ j=0  J ’ -I 

To  find  the  coefficients  a by  minimizing  H,  the  partial  derivative  of  H 
with  respect  to  each  a must  be  calculated,  the  resulting  expressions  set 
equal  to  zero,  and  the^a  's  solved  for  each  p=0,  1,  . . . n.  In  matrix 
form  the  solution  is:  ^ , 

A = B ■ -^R 


where: 


A .(ap) 
R "(r^) 
B 


p "*0 tl«  • • «n^ 


(Eq  25) 


k + L „ . . 

bj=  Z (t.)P^J. 
P'^  i = k - L ' 


3=0.  1. 


d.  Cautions  and  Conments.  The  equations  of  the  previous  paragraph 
niust  be  applied  with  care  when  prograimed  on  a computer.  If  the  size  of 
t^  is  large  relative  to  t^  - t^  _ p then  it  is  advantageous  to  translate 

the  t^'s  so  that  t^.  = 0 in  order  to  improve  numerical  accuracy.  After  a 

translation  it  can  be  observed  that  if  the  t^.  were  equally  spaced,  then 

B (already  a symmetric  matrix)  is  a constant  matrix,  a fact  that  could  be 
used  to  shorten  the  number  of  computations.  If  an  orthogonal  set  of  poly- 
nomials could  be  used  (only  reasonably  possible  in  the  case  of  equal 
spacing)  the  matrix  B would  be  diagonal,  and  formulas  of  the  type  mentioned 
in  subparagraph  a above  would  be  obtained.  Table  4 contains  the  result 
of  a trial  with  n=3,  m=15.  Figure  4 shows  a plot  of  data  joined  by 
straight  lines  and  a plot  of  smoothed  points  joined  by  straight  lines 
in  the  x-y  plane.  Figure  5 is  a plot  of  resultant  velocity  and  resultant 
acceleration  from  the  smoothed  values.  Further  information  is  contained 
in  Ralston. 
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Table  4.  Least  squares  fit  usinq  15  points,  degree  3 


X 

meters 

y 

meters)  ( 

X y z distance 

t(sec)  (meters)  (meters)  (meters 


37009.50  55188.278  79249.262  383.526  0 

37009.60  55188.695  79248.961  383.475  .5181 

37009.78  55189.771  79248.366  383.363  1.752 

37009.83  55190.086  79248.197  383.332  2.112 

37009.91  55190.594  79247.934  383.284  2.685 


distance  distance  between 
traveled  observed  & smoothed 


0 .1903 

.5181  .1149 

1.752  .1420 

2.112  .0351 

2.685  .0381 


t(sec) 


37009.50 
37009.60 
37009.78 
37009.83 
37009.91 


37009.50 

37009.60 

37009.78 

37009.83 

37009.91 


m/sec) 


y 

(m/sec) 


z 

(m/sec) 


vel oci ty 
resul  tant 


901 

- 3.071 

868 

- 3.163 

,021 

- 3.341 

233 

- 3.435 

,443 

-3.438 

;ec2) 

®y  2 
(m/sec  ) 

( 

(m, 

).754 

- 0.864 

- 0.994 

- 0.738 

- 0.612 

1 

L.625 

- 1 

1.644 

- .1 

L.385 

- 1 

L.083 

- 0.696 

- 1 

.588 

6.679 

.598 

6.694 

.604 

6.912 

.580 

7.141 

.644 

7.332 

c') 

acceleration 

resultant 

.147 

.909 

.803 

1.516 

1.288 
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6.  CONVOLUTION  FILTERS. 


a.  Transfoms.  The  objective  in  convolution  filtering  is  to  suppress 
the  changes  in  the  data  caused  by  errors  in  the  measuring  system  and  not 
to  change  the  desired  information.  In  the  ATMT  study  the  position  changed 
slowly  in  comparison  to  the  errors.  In  this  filtering,  the  original  data 
are  transformed  by  the  z-transform  from  a function  of  time  to  a function 
of  frequency.  If  T is  the  sampling  interval  (constant),  N the  number  of 

samples,  n = 2tt/NT,  and  j = / _ i , then  the  transformed  data  are 
X(kn): 

N - 1 

X(kfi)  = Z x(nT)e'J‘^^"'‘,  k=0,  1,  ...  N - 1. 
n=0 

Next,  a function  H(kn)  is  determined  such  that  H(kn)  is  close  to  1 for 
frequencies  Zirk/NT  near  the  frequencies  of  the  position  and  H(kP)  is  close 
to  0 for  frequencies  Zirk/NT  near  the  error  frequencies.  The  selection  of 
H(k52),  the  response  function,  is  discussed  in  the  next  subparagraph.  The 
functions  X(kn)  and  H(kfi)  are  multiplied  so  that  the  resulting  function 
is  approximately  X(kfi)  for  the  frequencies  of  the  position  and  approximately 
zero  for  the  frequencies  of  the  errors.  To  complete  the  filtering,  the 
inverse  z-transform  of  H*X  is  computed. 

N - 1 

y(nT)  = z X(kn)  H(kn)e 

^ R=0 


Thus,  the  main  problem  is  to  determine  N and  H(kn)  to  ensure  y(nT)  is  close 
to  the  true  position  part  of  x(nT). 

b.  Response  Functions.  A large  class  of  filters  may  be  derived  from 
HCz)  of  the  form: 


I 


H(z)  = i 

w (1  - e ^i  z ) 
i=l 


These  filters  include  Butterworth,  Chebyshev,  and  Bessel  filters  for 

different  choices  of  the  s^.  The  s^  for  the  Butterworth  filter  are  located 

equally  spaced  on  circle  in  the  s complex  plane  of  radius  w , w being 

c c 
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the  cutoff  frequency  of  lH(w | = .707.  The  s.  for  Chebyshev  filter  are 

w 1 

located  on  an  ellipse  in  the  s-plane.  Graphs  of  typical  response  functions 
H(2)  are  shown  in  fiqures  6 and  7. 


1 


Figure  6.  Butterworth  Response  Function  H(z)  or  |F(jw)|^  = 

1 + (w/w^)2n 


Figure  7.  Chebyshev  Response  Function  H(2)  or  |F(jw)|^  ^ 

7 — 5 

a Chebyshev  polynomial  of  order  n 1 + e (w/w^.) 
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Further  development  can  be  found  in  Gold  and  Rader,  reference  5.  Program 
versions  of  a third  order  Butterworth  filter  can  be  found  in  appendix 
(3  of  /\MSAA's  TR151,  reference  1. 


7.  MINIMUM  MAXIMUM  TECHNIQUES  FOR  FUNCTIONAL  APPROXIMATION.  This 
approach  is  generally  applied  when  the  value  of  a function  needs  to  be 
approximated  at  a large  number  of  domain  points  in  some  interval  and 
it  is  not  known  a priori  what  these  points  will  be.  Any  type  of 
approximating  functions  (e.g.,  polynomials,  rational  functions,  spline 
functions,  exponential  functions,  other  transcendental  functions)  can  be 
used,  but  rational  functions  are  the  type  used  most  often.  The  type  is 
determined  from  a judgment  of  which  type  is  most  "like"  the  data  being 
fitted  as  well  as  the  computational  ease  and  time  considerations.  Given 
the  type  of  function  to  be  used,  the  problem  in  this  method  is  to  choose 
the  function  of  the  given  type  that  minimizes  the  maximum  deviation 
between  the  observed  values  and  the  approximating  function.  For  the 
ATMT  application  it  was  unreasonable  to  expect  any  one  function  type 
(except  for  the  spline  function  approach)  to  be  able  to  follow  well  a 
vehicle's  path  for  2 to  3 minutes  without  using  large  degree  rational 
functions.  Working  the  problem  in  subintervals  or  being  concerned  with 
evaluation  at  all  points  is  similar  to  the  problems  discussed  in  paragraphs 
4 and  5.  For  further  reference,  see  Ralston,  reference  8,  or  Rice, 
reference  9. 

8.  SUMMARY.  The  decision  on  which  method  is  most  appropriate  for  use 
in  a given  application  must  be  based  on  judgment,  weighing  the  consider- 
ations discussed  in  the  previous  paragraphs.  Standard  rules  cannot  be 
developed  for  these  decisions;  however,  table  5 summarizes  some  of  the 
characteristics  of  the  various  approaches  that  should  be  considered. 
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Table  6.  Sunriary  of  method  characteristics 


1 

Better  for 

Sparse  (few) 
or  Dense  (many) 
Independent 
Variables 

Computation 

Ease 

Require 

Equal 

Spacing 

Derivatives 
Readily 
Available  at 
Points  other 
than  Data  Points 

1 Difference 
Quotients 

Neither 

Easy 

No 

No 

Interpolation 

Sparse 

Hard 

No 

Yes 

Kalman  Filter 

Dense 

Hard 

No 

No 

' Spline 

Sparse 

Hard 

No 

Yes 

Least  Squares 

! a.  Equal 

1 Spaces 

Either 

Easy 

Yes 

No 

1 b.  Not 

1 Equal 

‘ Spaces 

Either 

Hard 

No 

No 

Convolution 
' Filters 

Either 

Hard 

No 

No 

Min  Max 

Errors 

i 

k 

Sparse 



Hard 

No 

Depends 
on  Function 
Used 
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